Sixth Week: Linear Models

با توجه به داده های قیمت منازل لطفا با سوالات زیر پاسخ دهید.

۱. ماتریس همبستگی متغیرهای مختلف را به دست آورده و سپس رسم نمایید. اعداد به دست آمده را با آزمون فرض معناداری همبستگی بسنجید و سپس ده متغیری که همبستگی بالاتری با قیمت دارند را مشخص نمایید.

train %>% select_if(is.numeric)  -> train_nuemric
cormat <- round(cor(train_nuemric), 4)

melted_cormat <- melt(cormat, na.rm = TRUE)

ggplot(data = melted_cormat, aes(x = Var1, y = Var2, fill = value)) +
  geom_tile(color = "white") +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", 
                       midpoint = 0, limit = c(-1,1), space = "Lab",name="Pearson\nCorrelation")+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, 
                                   size = 8, hjust = 1),
        axis.text.y = element_text(angle = 0, vjust = 1, 
                                   size = 8, hjust = 1))+
  coord_fixed()

library("Hmisc")
res2 <- rcorr(as.matrix(train_nuemric))
# Extract the correlation coefficients
kable(res2$r, digits = 3)

# Extract p-values
kable(res2$P, digits = 3)

melted_cormat %>%
  filter(Var1 == "SalePrice") %>%
  filter(Var2 != "SalePrice") %>% 
  top_n(n = 10, wt = value) %>%
  arrange(-value) -> top_cor
kable(top_cor)

train_nuemric %>% select(SalePrice, OverallQual, GrLivArea, GarageCars, GarageArea,
                         TotalBsmtSF, `1stFlrSF`, FullBath, TotRmsAbvGrd, YearBuilt,
                         YearRemodAdd) -> train_top

۲. در یک تصویر نمودار پراکنش دو به دو ده متغیر بدست آمده به همراه قیمت را رسم نمایید و هم خطی بودن متغیرها را بررسی کنید

طبق نمودار ها می توان دید که برخی جفت ها با هم خطی دارند. مثل SalePrice-OverallQual, SalePrice-GrLivArea, GrLivArea-1stFlrSF, TotalBsmtSF-1stFlrSF, GrLivArea-TotRmsAbvGrd.

برای بررسی وجود هم خطی میان متغیرها از تست Farrar-Glauber و چند تست دیگر استفاده می کنیم. مقادیر این تست به ما نشان می دهند که برای مدل خطیمان هم خطی داریم که چیز خوبی نیست. نتیجتن باید چند تا از ستون ها را حذف کنیم تا هم خطی را از بین ببریم.

library(car)
scatterplotMatrix(train_top)

library(GGally)
ggpairs(train_top)

library(mctest)
omcdiag(train_top %>% select(-SalePrice), train_top$SalePrice)

## 
## Call:
## omcdiag(x = train_top %>% select(-SalePrice), y = train_top$SalePrice)
## 
## 
## Overall Multicollinearity Diagnostics
## 
##                        MC Results detection
## Determinant |X'X|:         0.0007         1
## Farrar Chi-Square:     10652.4387         1
## Red Indicator:             0.4852         0
## Sum of Lambda Inverse:    34.9056         0
## Theil's Method:           -0.2152         0
## Condition Number:        545.2532         1
## 
## 1 --> COLLINEARITY is detected by the test 
## 0 --> COLLINEARITY is not detected by the test

imcdiag(train_top %>% select(-SalePrice), train_top$SalePrice)

## 
## Call:
## imcdiag(x = train_top %>% select(-SalePrice), y = train_top$SalePrice)
## 
## 
## All Individual Multicollinearity Diagnostics Result
## 
##                 VIF    TOL       Wi       Fi Leamer    CVIF Klein
## OverallQual  2.7486 0.3638 281.7235 317.1576 0.6032 -0.2222     0
## GrLivArea    5.0188 0.1993 647.4656 728.9011 0.4464 -0.4056     1
## GarageCars   5.2510 0.1904 684.8768 771.0177 0.4364 -0.4244     1
## GarageArea   4.9310 0.2028 633.3298 712.9874 0.4503 -0.3986     1
## TotalBsmtSF  3.6027 0.2776 419.3196 472.0599 0.5269 -0.2912     0
## 1stFlrSF     3.6847 0.2714 432.5392 486.9422 0.5210 -0.2978     0
## FullBath     2.2155 0.4514 195.8305 220.4613 0.6718 -0.1791     0
## TotRmsAbvGrd 3.3565 0.2979 379.6611 427.4133 0.5458 -0.2713     0
## YearBuilt    2.3458 0.4263 216.8276 244.0993 0.6529 -0.1896     0
## YearRemodAdd 1.7510 0.5711 121.0020 136.2211 0.7557 -0.1415     0
## 
## 1 --> COLLINEARITY is detected by the test 
## 0 --> COLLINEARITY is not detected by the test
## 
## GarageArea , TotRmsAbvGrd , coefficient(s) are non-significant may be due to multicollinearity
## 
## R-square of y on all x: 0.7737 
## 
## * use method argument to check which regressors may be the reason of collinearity
## ===================================

۳. یک مدل خطی بر اساس ده متغیر برای پیش بینی قیمت برازش دهید. و سپس خلاصه نتایج مدل را به دست آورید.

fit <- lm(SalePrice ~ ., data=train_top)
summary(fit) # show results

## 
## Call:
## lm(formula = SalePrice ~ ., data = train_top)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -489958  -19316   -1948   16020  290558 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.186e+06  1.291e+05  -9.187  < 2e-16 ***
## OverallQual   1.960e+04  1.190e+03  16.472  < 2e-16 ***
## GrLivArea     5.130e+01  4.233e+00  12.119  < 2e-16 ***
## GarageCars    1.042e+04  3.044e+03   3.422 0.000639 ***
## GarageArea    1.495e+01  1.031e+01   1.450 0.147384    
## TotalBsmtSF   1.986e+01  4.295e+00   4.625 4.09e-06 ***
## `1stFlrSF`    1.417e+01  4.930e+00   2.875 0.004097 ** 
## FullBath     -6.791e+03  2.682e+03  -2.532 0.011457 *  
## TotRmsAbvGrd  3.310e+01  1.119e+03   0.030 0.976404    
## YearBuilt     2.682e+02  5.035e+01   5.328 1.15e-07 ***
## YearRemodAdd  2.965e+02  6.363e+01   4.659 3.47e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 37920 on 1449 degrees of freedom
## Multiple R-squared:  0.7737, Adjusted R-squared:  0.7721 
## F-statistic: 495.4 on 10 and 1449 DF,  p-value: < 2.2e-16

۴. نمودار قیمت واقعی و قیمت پیش بینی را رسم نمایید و خوب بودن مدل را ارزیابی کنید.

برای اینکه مدل خوب باشد باید داده ها نزدیک خط \(x=y\) باشند. که تا حد قابل قبولی این طور هست.

train_top %>% 
  mutate(PricePredict = predict(fit)) -> train_top_pre

ggplot(train_top_pre) +
  geom_point(aes(x = SalePrice, y = PricePredict), alpha = 0.3, colour = "blue")+
  geom_abline(slope = 1, intercept = 0, color = "red", size = 1, alpha = 0.5)

۵. مقدار R-squared مدل را به دست آورید. آیا بر اساس این کمیت مدل به خوبی به داده ها برازش داده شده است؟ کمیت F-statistic را در خلاصه مدل تفسیر نمایید.

مقدار R-squared مشخص می کند که مدلمان تا چه حدی تغییرات داده را می تواند توضیح دهد. مقدار adjusted-R-sqaured هم تاثیر زیاد بودن تعداد متغیرهای رگرسیون را کمتر می کند. این مقدار حدودن می شود \(0.77\) که تا حد خوبی قابل قبول است.

کمیت f-statistic اگر زیاد باشد نشانه ی وجود ارتباطی میان predictor های ما و response است. که در اینجا از آنجایی که مقدار آن زیاد است و پی ولیو هم بسیار کوچک است می توان این نتیجه را گرفت که ارتباطی وجود دارد. در واقع بیشتر بودن این مقدار اعتماد ما را به r-squared به دست آمده بیشتر می کند.

summary_fit<- summary(fit)
summary_fit

## 
## Call:
## lm(formula = SalePrice ~ ., data = train_top)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -489958  -19316   -1948   16020  290558 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.186e+06  1.291e+05  -9.187  < 2e-16 ***
## OverallQual   1.960e+04  1.190e+03  16.472  < 2e-16 ***
## GrLivArea     5.130e+01  4.233e+00  12.119  < 2e-16 ***
## GarageCars    1.042e+04  3.044e+03   3.422 0.000639 ***
## GarageArea    1.495e+01  1.031e+01   1.450 0.147384    
## TotalBsmtSF   1.986e+01  4.295e+00   4.625 4.09e-06 ***
## `1stFlrSF`    1.417e+01  4.930e+00   2.875 0.004097 ** 
## FullBath     -6.791e+03  2.682e+03  -2.532 0.011457 *  
## TotRmsAbvGrd  3.310e+01  1.119e+03   0.030 0.976404    
## YearBuilt     2.682e+02  5.035e+01   5.328 1.15e-07 ***
## YearRemodAdd  2.965e+02  6.363e+01   4.659 3.47e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 37920 on 1449 degrees of freedom
## Multiple R-squared:  0.7737, Adjusted R-squared:  0.7721 
## F-statistic: 495.4 on 10 and 1449 DF,  p-value: < 2.2e-16

۶. بر اساس p-value سطح معناداری ضرایب تصمیم بگیرید که چه متغیرهایی در مدل سازی استفاده شود. بر اساس متغیرهای جدید دوباره مدل سازی کنید و نتایج رو گزارش دهید.

بر اساس پی ولیوها به این نتیجه می رسیم که دو تا از عوامل ربط معناداری ندارند و آن ها را حذف می کنیم. حال اگر دوباره یک مدل خطی به داده هایمان فیت کنیم خواهیم دید که r-sqaured تغییر چندانی نمی کند یعنی مدل ما همان سطح از تغییرات را می تواند توجیه کند، ولی f-statistic ما بیشتر شده که بیانگر این است که ارتباط بیشتری بین predicatorها و response وجود است واطمینان ما از مقدار r-squared را بیشتر می کند.

train_top %>% 
  select(-GarageArea, -TotRmsAbvGrd) -> train_sig

fit_sig <- lm(SalePrice ~ ., data=train_sig)
summary(fit_sig)

## 
## Call:
## lm(formula = SalePrice ~ ., data = train_sig)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -482525  -19191   -1801   16208  289639 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.188e+06  1.284e+05  -9.255  < 2e-16 ***
## OverallQual   1.959e+04  1.188e+03  16.486  < 2e-16 ***
## GrLivArea     5.177e+01  3.097e+00  16.714  < 2e-16 ***
## GarageCars    1.395e+04  1.817e+03   7.680 2.92e-14 ***
## TotalBsmtSF   2.039e+01  4.269e+00   4.775 1.98e-06 ***
## `1stFlrSF`    1.465e+01  4.919e+00   2.979  0.00294 ** 
## FullBath     -7.184e+03  2.644e+03  -2.717  0.00666 ** 
## YearBuilt     2.699e+02  5.017e+01   5.380 8.69e-08 ***
## YearRemodAdd  2.957e+02  6.362e+01   4.649 3.64e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 37920 on 1451 degrees of freedom
## Multiple R-squared:  0.7734, Adjusted R-squared:  0.7721 
## F-statistic: 618.9 on 8 and 1451 DF,  p-value: < 2.2e-16

۷. مدل خود را بر اساس باقی مانده نقص یابی کنید. سه محک normality, independance, Constant Variance را در نقص یابی خود در نظر بگیرید.

Constant-Variance: برای بررسی این نکته از تست Breusch-Pagan استفاده می کنیم که فرض صفر آن این است که واریانس باقیمانده ها ثابت است. که رد می شود. پس این نکته درست نیست. همچنین با مشاهده ی نمودار باقی مانده ها و مقادیر فیت شده هم می توانیم ببینیم که یک روند وجود دارد و ثابت نیست.

independence: برای بررسی این نکته بررسی می کنیم تا ببینم آیا autocorrelation بین باقی مانده ها وجود دارد یا نه. برای بررسی این نکته از نمودار acf استفاده می کنیم. اگر باقی مانده ها autocorrelated نباشند باید بعد از lag0 نمودارمون بین ۲ خط آبی بیفتد (significance-level) که همین اتفاق هم افتاده. پس autocorrelation نداریم. راه دیگر برای بررسی این موضوع هم تست های randomness و durbin-watson اند. که نشان می دهند شواهد کافی برای رد کردن فرض autocorrelated نبودن وجود ندارد.

normality: برای بررسی این نکته از نمودار qqnorm استفاده می کنیم. اگر باقی مانده های ما روی یک خط باشند معنای آن این است که باقی مانده های ما نرمال هستند. راه دیگر برای بررسی این موضوع هم استفاده از تست های بررسی نرمال بودن است. در ادامه هم مشاهده می کنیم که طبق ۴ تا از این تست ها فرض نرمال بودن باقی مانده ها رد می شود. پس باقی مانده هایمان نرمال نیستند.

در آخر هم از تابع gvlma استفاده می کنیم که به طور خلاصه فرضیات رگرسیون خطی را برای مدلمان بررسی می کند.

پس در مجموع داده های ما در فرض های رگرسیون خطی صدق نمی کنند.

# Evaluate homoscedasticity
# non-constant error variance test
lmtest::bptest(fit_sig)  # Breusch-Pagan test

## 
##  studentized Breusch-Pagan test
## 
## data:  fit_sig
## BP = 395.73, df = 8, p-value < 2.2e-16

ncvTest(fit_sig)  # Breusch-Pagan test

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 2921.508    Df = 1     p = 0

# plot studentized residuals vs. fitted values 
spreadLevelPlot(fit_sig)

## 
## Suggested power transformation:  0.6280349

olsrr::ols_plot_resid_fit(fit_sig)

#No autocorrelation of residuals
acf(fit_sig$residuals)

# Method 2: Runs test to test for randomness
lawstat::runs.test(fit_sig$residuals)

## 
##  Runs Test - Two sided
## 
## data:  fit_sig$residuals
## Standardized Runs Statistic = -0.47124, p-value = 0.6375

# Method 3: Durbin-Watson test
lmtest::dwtest(fit_sig)

## 
##  Durbin-Watson test
## 
## data:  fit_sig
## DW = 1.9835, p-value = 0.3763
## alternative hypothesis: true autocorrelation is greater than 0

#=>   Durbin-Watson test


# Normality of Residuals
# qq plot for studentized resid
olsrr::ols_plot_resid_qq(fit_sig)

olsrr::ols_test_normality(fit_sig)

## -----------------------------------------------
##        Test             Statistic       pvalue  
## -----------------------------------------------
## Shapiro-Wilk              0.8062         0.0000 
## Kolmogorov-Smirnov        0.1059         0.0000 
## Cramer-von Mises         123.5906        0.0123 
## Anderson-Darling          35.91          0.0000 
## -----------------------------------------------

# all together
gvlma::gvlma(fit_sig)

## 
## Call:
## lm(formula = SalePrice ~ ., data = train_sig)
## 
## Coefficients:
##  (Intercept)   OverallQual     GrLivArea    GarageCars   TotalBsmtSF  
##   -1.188e+06     1.959e+04     5.177e+01     1.395e+04     2.038e+01  
##   `1stFlrSF`      FullBath     YearBuilt  YearRemodAdd  
##    1.465e+01    -7.184e+03     2.699e+02     2.957e+02

۸. داده ها را به پنج قسمت تقسیم کنید. بر اساس چهار قسمت مدل خطی را بسازید و صحت مدل را برای یک قسمت باقی مانده را تست کنید. خطای پیش بینی شما چقدر است؟

برای محاسبه ی خطا میانگین مربعات خطای حاصل از اختلاف مقدار پیش بینی شده و مقدار واقعی را به دست می آوریم. همچنین میانگین مربعی خطا ها (MSE) را هم به دست می آوریم.

set.seed(200)
#shuffling
train_sig_cv <- train_sig[sample(nrow(train_sig)),]

train_sig_cv_train <-
  train_sig_cv  %>% filter(row_number() < 0.8 * n())
train_sig_cv_test <-
  train_sig_cv  %>%  filter(row_number() >= 0.8 * n())

testmodel = lm(SalePrice ~ ., data= train_sig_cv_train)


train_sig_cv_test_pre_err <-
  train_sig_cv_test %>%
  mutate(pred =  predict.lm(testmodel, train_sig_cv_test, type = "response")) %>% 
  mutate(err = (SalePrice - pred) * (SalePrice - pred)) %>% 
  dplyr::select(SalePrice, pred, err)

model_error <- mean(train_sig_cv_test_pre_err$err)
model_error

## [1] 1148541154

sqrt(model_error)

## [1] 33890.13

۹. آیا قیمت ربط غیر خطی با یکی از ده متغیر استفاده شده دارد؟ بر اساس دستاوردهای خود مدل را بهتر نمایید.

نمودار SalePrice را در مقابل سایر متغیرها می کشیم و مشاهده می کنیم که ۲ تا از متغیرها ربط غیرخطی دارند و به نظر از درجه ی ۲ دارند.

ولی برای بهترشدن مدل و کسب رتبه ی بهتر در kaggle! به روش آزمون و خطا از درجات بالاتر سایر متغیرها هم استفاده می کنیم تا بهتر بتوانیم آن ها را توصیف کنیم. همچنین به جای خود SalePrice از لگاریتم آن استفاده می کنیم.

sp <- train_top$SalePrice
for(i in 2:11){
  other <- train_top[, i]
  sp_other <- data.frame(sp = sp,
                   xaxis = other)
  colnames(sp_other) <- c("sp", "xaxis")
  ggplot(sp_other, aes(x=xaxis, y=sp)) + 
    geom_point(alpha = 0.3)+
    geom_smooth(method = "loess")+
    geom_smooth(method = "lm", color = "red", alpha = 0.5)+
    xlab(colnames(other))+
    ylab(colnames(sp)) -> graph
  print(graph)
}

#####
#####   OverallQual
other <- train_top[, 2]
sp_other <- data.frame(sp = sp,
                 xaxis = other*other)
colnames(sp_other) <- c("sp", "xaxis")
ggplot(sp_other, aes(x=xaxis, y=sp)) + 
  geom_point(alpha = 0.3)+
  geom_smooth(method = "loess")+
  geom_smooth(method = "lm", color = "red", alpha = 0.5)+
  xlab(paste(colnames(other), "^2"))+
  ylab(colnames(sp)) -> graph
print(graph)

#####
#####   TotalBsmtSF
other <- train_top[, 6]
sp_other <- data.frame(sp = sp,
                       xaxis = other *other)
colnames(sp_other) <- c("sp", "xaxis")
ggplot(sp_other %>% filter(xaxis < 30000000), aes(x=xaxis, y=sp)) + 
  geom_point(alpha = 0.3)+
  geom_smooth(method = "loess")+
  geom_smooth(method = "lm", color = "red", alpha = 0.5)+
  xlab(paste(colnames(other), "^2"))+
  ylab(colnames(sp)) -> graph
print(graph)

#####
#####   improving model

train_top_improved <- train_top %>% 
  mutate(TotalBsmtSF2 = TotalBsmtSF ** 2,
         TotalBsmtSF3 = TotalBsmtSF ** 3,
         OverallQual2 = OverallQual ** 2,
         OverallQual3 = OverallQual ** 3,
         GrLivArea2 = GrLivArea ** 2,
         GarageArea2 = GarageArea ** 2,
         YearBuilt2 = YearBuilt ** 2)

model_improved <- lm(log2(SalePrice) ~ ., data = train_top_improved %>%
                       dplyr::select(-TotRmsAbvGrd, -GarageArea, -TotalBsmtSF, -`1stFlrSF`))
summary(model_improved)

## 
## Call:
## lm(formula = log2(SalePrice) ~ ., data = train_top_improved %>% 
##     dplyr::select(-TotRmsAbvGrd, -GarageArea, -TotalBsmtSF, -`1stFlrSF`))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.84406 -0.10406  0.01634  0.12898  0.80179 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.748e+02  2.947e+01  -5.932 3.73e-09 ***
## OverallQual   4.367e-01  9.351e-02   4.670 3.29e-06 ***
## GrLivArea     6.943e-04  5.298e-05  13.105  < 2e-16 ***
## GarageCars    8.263e-02  1.411e-02   5.858 5.78e-09 ***
## FullBath     -4.828e-02  1.608e-02  -3.003 0.002715 ** 
## YearBuilt     1.826e-01  2.990e-02   6.108 1.30e-09 ***
## YearRemodAdd  4.258e-03  4.047e-04  10.521  < 2e-16 ***
## TotalBsmtSF2  1.681e-07  1.183e-08  14.216  < 2e-16 ***
## TotalBsmtSF3 -3.499e-11  2.569e-12 -13.619  < 2e-16 ***
## OverallQual2 -5.279e-02  1.575e-02  -3.352 0.000822 ***
## OverallQual3  2.908e-03  8.565e-04   3.396 0.000703 ***
## GrLivArea2   -7.637e-08  1.338e-08  -5.707 1.40e-08 ***
## GarageArea2   1.066e-07  4.475e-08   2.381 0.017373 *  
## YearBuilt2   -4.598e-05  7.658e-06  -6.005 2.42e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2226 on 1446 degrees of freedom
## Multiple R-squared:  0.8521, Adjusted R-squared:  0.8508 
## F-statistic: 640.8 on 13 and 1446 DF,  p-value: < 2.2e-16

۱۰. بر اساس مدل نهایی به دست آمده نتایج پیش بینی خود را بر روی test.csv به دست آورید و در سایت kaggle در مسابقه House Prices: Advanced Regression Techniques بارگذاری نمایید. سپس لینک رتبه و عدد آن را ضمیمه تمرین کنید.

read_csv("./../house/test.csv") -> test_data

test_data_2 <- test_data %>% 
  mutate(TotalBsmtSF2 = TotalBsmtSF ** 2,
         TotalBsmtSF3 = TotalBsmtSF ** 3,
         OverallQual2 = OverallQual ** 2,
         OverallQual3 = OverallQual ** 3,
         OverallQual4 = OverallQual ** 3,
         GrLivArea2 = GrLivArea ** 2,
         GarageArea2 = GarageArea ** 2,
         YearBuilt2 = YearBuilt ** 2)
test_data_2[is.na(test_data_2)] <- 0

test_data_2 %>% 
  mutate(SalePrice = 2 ** predict.lm(model_improved, test_data_2)) %>% 
  dplyr::select(Id, SalePrice) -> result

write_csv(result, "./../result.csv")

نتیجه و رتبه

score